A proof-of-concept study on the genomic evolution of Sars-Cov-2 in molnupiravir-treated, paxlovid-treated and drug-naïve patients

Little is known about SARS-CoV-2 evolution under Molnupiravir and Paxlovid, the only antivirals approved for COVID-19 treatment. By investigating SARS-CoV-2 variability in 8 Molnupiravir-treated, 7 Paxlovid-treated and 5 drug-naïve individuals at 4 time-points (Days 0-2-5-7), a higher genetic distance is found under Molnupiravir pressure compared to Paxlovid and no-drug pressure (nucleotide-substitutions/site mean±Standard error: 18.7 × 10−4 ± 2.1 × 10−4 vs. 3.3 × 10−4 ± 0.8 × 10−4 vs. 3.1 × 10−4 ± 0.8 × 10−4, P = 0.0003), peaking between Day 2 and 5. Molnupiravir drives the emergence of more G-A and C-T transitions than other mutations (P = 0.031). SARS-CoV-2 selective evolution under Molnupiravir pressure does not differ from that under Paxlovid or no-drug pressure, except for orf8 (dN > dS, P = 0.001); few amino acid mutations are enriched at specific sites. No RNA-dependent RNA polymerase (RdRp) or main proteases (Mpro) mutations conferring resistance to Molnupiravir or Paxlovid are found. This proof-of-concept study defines the SARS-CoV-2 within-host evolution during antiviral treatment, confirming higher in vivo variability induced by Molnupiravir compared to Paxlovid and drug-naive, albeit not resulting in apparent mutation selection.

Molnupiravir is a bioactive isopropylester prodrug of the ribonucleoside analogue β-D-N4-hydroxycytidine (NHC, EIDD1931) used as a broad-spectrum antiviral agent and originally developed against alphaviruses. Following administration, Molnupiravir is converted into active NHC in the plasma by host esterases and is then intracellularly phosphorylated by host kinases, leading to its triphosphate form (NHC-TP) 1 . In this form, it can act as an alternate and competitive substrate for the SARS-CoV-2 RNA-dependent RNA polymerase (RdRp), interfering with viral replication. Similar to other mutagenic nucleotides, NHC-TP can affect base pairing in various ways, due to the existence of different tautomeric forms: in its hydroxylamine form, NHC-TP acts like cytosine (C), thus leading to base pairing with guanine (G), while in its oxime form acts like uracil (U), leading to adenine (A) base pairing. During the sequential steps of RNA replication, the RdRp incorporates NHC-TP into the negative or positive sense genomic RNA (±gRNA) of SARS-CoV-2 as NHC-monophosphate (NHC-MP) 2,3 . When present as NHC-MP, both tautomeric forms coexist, resulting in the incorporation of GTP and ATP, causing an increase in G to A and C to U transitions in the final +gRNA template. As the viral RNA genome is amplified in the cell, the effect of this accumulation of mutations comprise mutagenic transcription and translation products, ultimately yielding nonfunctional genomes, in a process which is known as lethal mutagenesis or error catastrophe 4 . Differently from what happens for the nucleoside analogue Remdesivir, whose translocation barrier can in part be eliminated by the proofreading activity of the Nsp14 exonuclease 5 , the mutations induced by Molnupiravir are not corrected, due to the stability of the M:G and M:A base pairs formed by the NHC-MP 6 , thus resulting in higher antiviral activity.
Paxlovid (PF-07321332) is a co-packaged combination of Nirmatrelvir, a novel main protease (Mpro) inhibitor, specifically designed to block the activity of the SARS-CoV-2 Mpro, and Ritonavir, a strong cytochrome P450 (CYP) 3A4 inhibitor 7,8 . Paxlovid works intracellularly by binding to the highly conserved main protease of the SARS-CoV-2. By this mechanism, Paxlovid inhibits viral replication at the polyprotein maturation step, which occurs before viral RNA replication.
Data resulting from clinical trials have confirmed both Molnupiravir and Paxlovid antiviral efficacy and tolerability, highlighting their efficacy in reducing hospitalization and death in non-hospitalized adults with mild-to-moderate COVID-19, and reporting the absence of any major adverse effects [9][10][11] .
From a virological point of view, a significant decrease in the time required for viral RNA clearance and a higher proportion of overall viral RNA clearance was observed in participants treated with Molnupiravir compared to those not treated 12 . Similar data are predicted for Paxlovid 13 .
So far, limited information is available about Molnupiravir and Paxlovid mutagenic potential in patients, as well as about the risk of drug resistance emergence. It has been proposed that Molnupiravir could have a high genetic barrier to resistance 14,15 , due to the random accumulation of mutations throughout the genome, together with its short-term use, which suggests a low likelihood for the appearance of resistance as recently confirmed by the AGILE phase IIa clinical trial 16 . Regarding Paxlovid, cell culture systems revealed the selection of Nirmatrelvir resistant Mpro variants characterized by high fitness and by weakness in Nirmatrelvir-Mpro binding [17][18][19] .
In such a rapidly evolving pandemic, real-world data are now demanded to understand the SARS-CoV-2 evolution under small antiviral agents' pressure. In this light, this proof-of-concept study defines the dynamics of the evolution of SARS-CoV-2 genome in infected individuals across the five days of small antiviral drug treatment (Molnupiravir or Paxlovid), focusing the attention on the emergence of single nucleotide polymorphisms (SNPs) in the whole genome, and their persistence after the end of treatment (Day 7). Results are compared with the SARS-CoV-2 evolution observed in drug-naïve individuals. Among the 3 study groups, age and prevalence of comorbidities were higher in Molnupiravir and Paxlovid treated individuals compared to drug-naïve (P = 0.005 and 0.032). Molnupiravirtreated and Paxlovid-treated patients had a first positivity date more recent compared to drug-naïve individuals (P = 0.033), and the prevalence of the sublineage BA.2.9 was more prevalent in Molnupiravir-treated and drug-naïve individuals compared to Paxlovid-treated patients (P = 0.002) ( Table 1).  10 copies/mL respectively, P = 0.056), highlighting a more rapid decay kinetics with nucleoside analogue than what seen in Paxlovid-treated and drug-naïve individuals (Fig. 1).

SARS-CoV
The viral load observed at each time point was summarized below. At the beginning of Molnupiravir and Paxlovid treatment (Day 0), median SARS-CoV-2 nasopharyngeal load was 7.2 (IQR: 6.5-7.7) log copies/mL, and no significant differences were observed among Molnupiravir-treated, Paxlovid-treated and drug-naïve individuals (Table 1 and Fig. 1).
Within-host SARS-CoV-2 genetic distance and correlation with treatment. SARS-CoV-2 sequences obtained at different time points from Molnupiravir-treated, Paxlovid-treated and drug-naive individuals were subjected to phylogenetic analysis in order to define the SARS-CoV-2 evolution across time points and against treatment (Fig. 2).
ML tree first revealed that all sequences belonging to the same patient correctly clustered together, without intermixing, excluding cross-contamination among samples (Fig. 2). ML tree also confirmed that all sequences belonged to Omicron clade, sublineage BA.2.
Orf8, known to be involved in immune response evasion 24 , seemed to be the only protein evolving under a positive selection during Molnupiravir treatment (dN>dS, P = 0.001), due to the presence of 7/7 (100%) non-synonymous mutations. Among the mutated amino acid positions, none was under purifying selection, and only orf8 A27927G-T12A persisted between Day 5 (intra-patient prevalence: 10.3%) and Day 7 (intra-patient prevalence: 55.8%) of ID10.
When considering RdRp and Mpro, target sites of Molnupiravir and Paxlovid, respectively, no evidence of evolution under positive selection was reported (Table 3). Among the 39 SNPs inducing amino acid changes in the RdRp across time points in Molnupiravir-treated patients, 12 were localized in the fingers, 9 in the interface domain, 5 in the palm, and 3 in the thumb domain (Fig. 4). The remaining mutations were localized in the NiRAN domain and in the hairpin like structure. Among these missense mutations, no evidence of persistence was reported, except for G13604A-R55H appearing at Day 2 (intra-patient prevalence: 5.3%) and persisting at Day 5 (intra-patient prevalence: 14.3%) of ID11. No amino acid change was observed in the RdRp active site residues 25 .
Of note, no evidence of extensive missense mutational rate was reported in Mpro during Paxlovid treatment as well as no evidence of so far known (https://www.fda.gov/media/155050/ download) 17 Mpro mutations conferring resistance to Paxlovid was found in Paxlovid treated patients. Only the following Mpro mutations were found: C10301T-Q83* found at Day 7 of patient 4 with a relative abundance of 44.6%, and T10504C-F150F found at Day 2 of patient 3 with a relative abundance of 15.1%.
Interestingly, the mutation C10376T-P108S, whose emergence was observed during Nirmatrelvir treatment in SARS-CoV-2 cell culture, was observed at Day 5 of Molnupiravir treatment in patient 11 (relative abundance of 14.5%). However, no data of the reduction of Nirmatrelvir activity was reported for this mutation in biochemical assays (https://www.fda.gov/media/155050/ download). Since the loss of the nsp14 could increase the sensitivity of SARS-CoV-2 to Molnupiravir, we further investigate the effect of Molnupiravir in inducing variability on the proofreading function of nsp14. No evidence of evolution under positive selection was reported (Table 3). Among the 23 missense mutations, only two increased their prevalence from 11.1% and 5.4% at Day 5 to 52.7% and 45.3% at Day 7 of ID10 (nsp14 C18452T-A138V and G18520A-V161I). No amino acid change was observed in the nsp14 active site residues 26 .

Discussion
This proof-of-concept study on the SARS-CoV-2 genomic evolution under small molecule antivirals pressure indicates that Molnupiravir induced a higher within-host SARS-CoV-2 variability compared to Paxlovid, finding supported by the higher number of SNPs observed across time points. The SNPs induced by Molnupiravir appeared randomly along the SARS-CoV-2 genome, were most G to A and C to T transitions and showed their peak at Day 5 of treatment, when the genetic distance showed its highest value, and the viral load reached its lowest peak. Two days after the end of treatment (Day 7) the estimates of genetic distance, number of SNPs and viral load were stable compared to those reported at the end of treatment, probably reflecting the fact that, even if the treatment was stopped at Day 5, the antiviral drug was still present as a pharmacological tail and continued to drive the genetic evolution. These findings are consistent with the known error catastrophe induced by Molnupiravir against the virus 4 and are in line with recent large-scale data. In the AGILE phase IIa clinical trial 16 , an increase in G to A, C to T and T to C mutations was observed during Molnupiravir treatment, and no amino acid substitutions were enriched consistently at specific sites during the first 5 days of treatment.
Despite the intra-host SARS-CoV-2 variation, already observed among clinical samples throughout the pandemic 27 CoV-2 strains under the Molnupiravir and Paxlovid pressure, as well as in drug-naïve samples, did not seem to experience strong selective evolution. By calculating the probability to evolve under purifying selection of each gene, only orf8, known to be involved in immune response evasion 24 , seemed to evolve under a positive selection during Molnupiravir treatment, as a consequence of the 7/7 non-synonymous mutations detected randomly along the protein. Among these, only the orf8 A27927G-T12A persisted between Day 5 and Day 7. No clear information is available about the role of this amino acid mutation in protein stability or functionality, but in silico prediction suggests that this mutation might induce a protein stability change (free energy change, ΔΔG, of −1.5 Kcal/mol compared to wild type) 29 . Further in vitro studies are needed to confirm these in silico findings as well as in general to firmly establish whether and how these substitutions evolve under selection pressure.
We noted that this amino acid change was detected with other 4 amino acid mutations (the nsp3 G5572A-M951I and T8004C-V1762A, nsp14 C18452T-A138V and G18520A-V161I) at Day 5 of the Molnupiravir-treated patient ID10, a renal transplant recipient on tacrolimus-based immunosuppression. This mutational pathway was selected at Day 5 and spread at Day 7 when an increase of genetic distance and a posttreatment increase of viral load were detected. Among the mutations observed, the nsp14 mutations fall in the exonuclease domain, and thus might impair the functionality of the errorcorrection mechanism, leading to accumulation of mutations and increase of variability. These hypotheses require further validation by in vitro studies and docking simulations. The nsp3 G5572A-M951I mutation falls in the PLpro domain, which plays an essential role in the viral replication and the innate immune evasion 30,31 . If the substitution of methionine to Fig. 3 Circos plots representing the SNPs found in SARS-CoV-2 strains. Panels represent SARS-CoV-2 strains found in (a) Molnupiravir-treated, (b) Paxlovid-treated and (c) drug-naïve individuals across time points. Each time point is defined by one track and grouped according to patient. Tracks in black represent timepoints not available; striped tracks represent samples not sequenced as viral load was undetectable. SNPs are defined by different colors, shapes and size according to type of nucleotide substitution (red: G-A; orange: C-T; green: other), non-synonymus (circle) or synonymus mutations (triangle), and intra-patient prevalence, respectively. The outside track corresponds to SARS-CoV-2 genome architecture. Overlapping orfs are not reported and comprise: orf2b (coordinates 21744-21860) overlapping S; orf3b (coordinates 25814-25879), orf3c (coordinates 25457-25579), orf3d (coordinates 25524-25694) and orf3d-2 (coordinates 25596-25694) overlapping orf3a; orf9b (coordinates 28284-28574) and orf9c (coordinates 28734-28952) overlapping N. Table 3 Codon-based Neutrality for each SARS-CoV-2 gene against study groups.

SARS-CoV-2 proteins
Molnupiravir Paxlovid Drug naive Codon-based Test of Neutrality for each SARS-CoV-2 protein against study group. The test statistic (dN-dS) is reported as V, where dS and dN are the numbers of synonymous and non-synonymous substitutions per site, respectively. The variance of the difference was computed using the bootstrap method (500 replicates). The probability of rejecting the null hypothesis of strict-neutrality (dN = dS) is shown as P-values (P). P-values less than 0.05 are considered significant at the 5% level and are highlighted. Analyses were conducted using the Nei-Gojobori method. All ambiguous positions were removed for each sequence pair (pairwise deletion option). Evolutionary analyses were conducted in MEGA11.
isoleucine in PLpro secondary structure evoked effect on the protein stability, needs to be investigated. A progressive SARS-CoV-2 genetic variability over time points was observed also in ID12, but with different features. At Day 2 of Molnupiravir-treatment, two nsp14 non-synonymous mutations appeared (C18511A-P158T and C18647T-P203L) with an intrapatient prevalence of 36.0% and 23.7%, respectively. These mutations are both localized in the ExoN domain of nsp14, and the P203L was recently suggested to induce a weakly interaction between nsp14 and nsp10, possibly resulting in reduction of the proofreading activity of nsp14 32 . As a result, nsp14-P203L variants showed higher nucleotide substitution rate at population level, and in our patient induced a pronounced genomic evolution, confirmed by genetic distance (65.0 [±3.2] × 10 −4 ) and emergence of 175 minority SNPs (relative abundance %, median IQR 8.7 [6.5-13.6]) at Day 5.
Understanding the factors underlining the emergence of these mutational pathways in these two patients and whether they are correlated to health status rather than drug pressure might provide essential insights into the risk of SARS-CoV-2 evolutionary mechanisms.
By our proof-of-concept study, no clear evidence of RdRp and Mpro mutations that might contribute to Molnupiravir or Paxlovid resistance (https://www.fda.gov/media/155050/download) was found, even if further insights are needed to define if some of the enriched mutations might be involved in drug resistance mechanisms.
Looking at the variability of the whole spike protein, no site on positive selection or persistent mutation were detected. Six non-synonymous SNPs (V327I, C391F, V401I, T430I, L455S, and R457K) with a relative abundance of median (IQR) 9.8 (6.4-14.9) appeared at different time points in the RBD, but no chance of further selection or spreading was revealed, suggesting no evident risk of SARS-CoV-2 spike evolution under Molnupiravir.
Worthy of mention is the potential impact of the Molnupiravir-induced variability on molecular and Ag diagnostic assays. As already known, SARS-CoV-2 variability in diagnostics targets could hamper the viral load detection or Ag positivity 33 . Regarding the antigenic tests, false negative results might be considered unlikely events under Molnupiravir treatment, because of the rare persistence of non-synonymous mutations and the absence of positive selection in the target genes of antigen tests. With the molecular tests, instead, the risk of detection failure could be due to SNPs falling in the target region of the probes or primers, thus interfering with their binding and thus rendering the assays ineffective, but the multi-target RT-PCR tests used in clinical practice withdraws this risk. These hypotheses require further analyses to be in depth assessed and defined.
Another aspect that should be taken into consideration is the use of Molnupiravir in combination with Paxlovid-or monoclonal antibodies-based therapy. Since the mode of action of these molecules is different, combined therapy might represent an advantage compared to monotherapy, allowing to boost the efficacy against SARS-CoV-2 infection and to reduce the risk of drug resistance or immune escape mutations. Studies reporting a synergistic effect of Molnupiravir and Paxlovid combination strategy both in vitro 34,35 and in animal models 36 are already published, but ad hoc designed clinical studies will be needed to confirm these findings in humans.
Our study has some limitations. The analysis of negative and positive selection induced by Molnupiravir and Paxlovid should be interpreted carefully, as the number of SNPs that define the purifying selection is small and potential errors introduced during reverse transcription, PCR amplification, or sequencing can slightly impact the results 37 . To overcome these problems, only SNPs having a minimum supporting read frequency of 2% with a depth ≥20, and thus expected to provide more robust results, were retained.
The sample size was relatively small, limiting general considerations. We cannot exclude that denser sampling would reveal novel patterns of mutations not observed here and further address genetic variability. The small sample size might have precluded the possibility to detect drug resistance, whose emergence under small drugs therapies remains a challenge to investigate. The results obtained were mainly addressed to hospitalized and immunocompromised individuals. This population bias may have also caused a substantial underestimation of the viral load decay kinetics and a pronounced genomic variability across time.
In conclusion, this study contributes to define the dynamics of within-host variability of SARS-CoV-2 during Molnupiravir and Paxlovid treatment, confirming the higher genomic variation induced by this nucleoside analogue (as its main mechanism of action also in vivo), and describing in detail the effect of the drug upon the evolution of each different gene constituting the SARS-CoV-2 genome. Molnupiravir induces a number of SNPs that increase against days of treatment. Due to bottleneck or purifying selection, these SNPs appear randomly and rarely persist across time points, limiting the potential risk of selection of viral variants characterized by a fitness advantage during Molnupiravirtreatment. Ad hoc designed studies are necessary to confirm our findings, provide an extensive overview of SARS-CoV-2 intrahost variability and minority variants description during small molecule antivirals treatment, if and how these minority variants can spread in the population, and their potential role in virulence and transmissibility. (i) mutations appearing at day 2 with relative abundance between 2 and 10% in light green, and between 10 and 20% in green; (ii) mutations appearing at day 5 with relative abundance between 2 and 10% in light blue, and between 10 and 20% in blue; (iii) mutations appearing at day 7 with relative abundance between 2 and 10% in light red, between 10 and 20% in red, and between 20 and 40% in dark red; (iv) mutations appearing at both day 2 and day 5 in cyan (including the G13604A-R55H mutation); (v) mutations appearing at day 2 or day 7 in yellow; vi) mutations appearing at day 5 or day 7 in purple. Template RNA is colored in black, while product RNA in olive. RdRp: RNA-dependent RNA polymerase.

Methods
Study population. This retrospective observational proof-of-concept study initially included 24 individuals diagnosed for SARS-CoV-2 infection at Azienda Ospedaliero-Universitaria Policlinico of Modena (Italy) and Bambino Gesù Children Hospital IRCCS of Rome (Italy) from March 2022 to May 2022. Of these, 10 individuals started SARS-CoV-2 antiviral treatment with Molnupiravir, 9 with Paxlovid, and 5 individuals received no treatment (drug-naïve), depending on physicians' discretion.
For each individual with nasopharyngeal swabs available, SARS-CoV-2 viral load and genomic variability were retrospectively analyzed at 4 time points: before SARS-CoV-2 antiviral treatment start (Day 0), 2 days after SARS-CoV-2 antiviral treatment start (Day 2), 5 days after SARS-CoV-2 antiviral treatment start, i.e. end of treatment (Day 5), and 2 days after treatment end (Day 7). The same data were obtained at 4 comparable time points for the 5 drug-naive individuals ( Supplementary Fig. 1).
Based on the availability of samples at each time point, results on SARS-CoV-2 viral load quantification and whole-genome SARS-CoV-2 sequencing, 20 individuals were finally selected for the study, and divided in three study groups according to treatment as follows: 8 Molnupiravir-treated, 7 Paxlovid-treated and 5 drug-naïve individuals.
Demographic, epidemiological, and clinical data were obtained retrospectively from pseudonymized electronic medical records.
Ethics. The study protocol was approved by local Research Ethics Committee of the Azienda Ospedaliero-Universitaria Policlinico of Modena and Bambino Gesù Children hospitals (prot. 396/2020/OSS/AOUMO and prot. 2384_OPBG_2021). This study was conducted in accordance with the principles of the 1964 Declaration of Helsinki. Informed consent was waived in accordance with the hospital regulations on retrospective observational studies.
SARS-CoV-2 viral load quantification and sequencing. Total RNA was extracted from 280 µl of nasopharyngeal swabs (COPAN FLOQSwab with 3 mL of Universal Transport Medium) using QIAamp viral RNA mini kit (Qiagen) following the manufacturer's instruction. SARS-CoV-2 RNA was quantified by means of the QX200™ Droplet Digital™ PCR System (ddPCR, Bio-Rad Laboratories, Inc.) using home-made protocol targeting RdRp and the RNAseP housekeeping gene 38 . Quantitative results were then normalized in number copies/mL of swab.
Whole-genome sequences of SARS-CoV-2 were generated with 50 ng viral RNA template, by using the QIAseq DIRECT SARS-CoV-2 Kit according to the manufacturer's protocol 39,40 . Libraries were then sequenced on a MiSeq instrument (Illumina, San Diego, CA, USA) with 2 × 150-bp paired-end reads. Raw reads were trimmed for adapters and filtered for quality (average q28 threshold and read length >135 nt) using Fastp (v0.20.1) 41 . First and last 15 nucleotides were also removed from all reads. Reference-based assembly was performed with BWA-mem algorithm 42 , aligning against the GenBank reference genome NC_045512.2 (Wuhan, collection date: December 2019). SNPs variants were called through a pipeline based on samtools/bcftools 43 , and all SNPs having a minimum supporting read frequency of 2% with a depth ≥20 were retained. Consensus was generated using the github freely distributed software vcf_consensus_builder (https://github. com/peterk87/vcf_consensus_builder).
Genetic distance and Phylogenetic analysis. SARS-CoV-2 lineages of the SARS-CoV-2 consensus sequences obtained were assigned according to Pangolin application (Pangolin v4.1.1, https://github.com/cov-lineages/pangolin/releases). Sequences were aligned using MAFFT v7.475 and manually inspected using Bioedit. The final alignment comprised 65 sequences of 29,184 nucleotides of length. The within-host viral genetic distance (always expressed as nucleotide substitutions per site) over all time points, and between time points intervals (Day 0-Day 2; Day 2-Day 5; Day 5-Day 7) was calculated with MEGA (v11) 44 under the Tamura-Nei model 45 . Differences in average SARS-CoV-2 genetic distance among the study's groups was then evaluated by Kruskal-Wallis test. Differences in viral genetic distance between time points intervals within each study group were assessed by Wilcoxon signed-rank test.
In order to further explore the within-host and within-group evolution of SARS-CoV-2, a maximum likelihood (ML) phylogeny tree was performed with IqTree2 (v2.1.3) 46 with 1000 bootstrap replicates, using the best-fit model of nucleotide substitution TN + F + G4 inferred by ModelFinder 47 . Annotation of the phylogenetic tree, including information about study group, lineages, SARS-CoV-2 viral load and timepoints, was performed with iTOL 48 .
Non-synonymous and synonymous substitution rate estimation. We further analyzed genetic variability across all the 33 SARS-CoV-2 genes for evidence of purifying selection or neutral selection. Genes and their overlapping open reading frames (ORFs) were considered according to the reference genome (GenBank accession number NC_045512) and Jungreis et al. 2021 49 . The rates of synonymous (dS) and non-synonymous (dN) mutations per SARS-CoV-2 gene were computed on MEGA (v11), using the Nei-Gojobori method (Proportion) 50 , by removing all the constitutive SNPs and using the bootstrap method with 500 replicates. If genes evolve under neutral selection, the SNPs are likely to distribute equally at each codon position and the difference between non-synonymous substitution and synonymous substitution tends to be not significant (P > 0.05). On the contrary, positive (dN>dS) or negative (dN<dS) selection suggests that the SNPs are not equally distributed at each codon position and the difference between nonsynonymous substitution and synonymous substitution tends to be significant (P ≤ 0.05). In order to identify single sites under purifying selection (P ≤ 0.1) the Fixed Effects Likelihood (FEL) was applied 51 .
Statistics and reproducibility. Descriptive statistics are expressed as median values and interquartile range (IQR) for continuous data and number (percentage) for categorical data. To assess significant differences, chi-squared test for trend and Kruskal-Wallis were used for categorical and continuous variables, respectively. Jonckheere-Terpstra test was used to define differences in viral load, viral load decay and SARS-CoV-2 genetic variability at different time points among study groups.
The presence of SNPs at the different time points, together with the information about the type of mutation (synonymous vs non-synonymous), the type of nucleotide substitution (G-A, C-T, or other), and their relative abundance, were graphically represented as circus plots using the Circos package v0.69-9 52 .
Statistical analyses were performed with SPSS software package for Windows (version 23.0, SPSS Inc., Chicago, IL). A two-sided p-value <0.05 was considered statistically significant.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.